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In the last few years, the Imaging Atmospheric Cherenkov Telescopes H.E.S.S., MAGIC 
and VERITAS have detected more than 40 blazars in the very-high-energy range (VHE, 

100 GeV—100 TeV). During their trip to us, the VHE photons undergo an energy-dependent 
absorption by scattering off the infrared/optical/ultraviolet photons of the extragalactic 
background light (EBL), which is produced by galaxies during the whole cosmic evolution. 
Actually, both the observed spectra and the emitted ones predicted by conventional VHE 
photon emission models have a simple power-law behavior to a good approximation. Within 
conventional physics, the emitted slope T em (;r) changes randomly with the source redshift 2 
by a large amount. But, surprisingly, the {r em (z)} distribution exhibits a correlation with 
z, since the associated best-fit regression line is a decreasing function of z, leading blazars 
with harder spectra to be found only at larger redshift. So, how can the source distribution 
get to know the redshifts in such a way to adjust their T em (z) values so as to reproduce 
such a statistical correlation? It is very difficult to imagine an intrinsic mechanism which 
could lead to this spectral variation within conventional physics, given the fact that neither 
cosmological evolutionary effects nor observational selection effects can explain it. Things 
are quite different in the presence of axion-like particles (ALPs), which are predicted by 
several extensions of the Standard Model and especially by those based on superstring the¬ 
ories, and are attracting growing interest being also good candidates for cold dark matter. 

We show that photon-ALP oscillations occurring in extragalactic magnetic fields yield, for 
a realistic choice of the parameters, a {r em (-^)} distribution whose best-fit regression line 
becomes amazingly redshift-independent - indeed in agreement with the physical intuition 
- and so the above problems disappear. This is evidently a highly nontrivial fact, which 
therefore provides preliminary evidence for the existence of ALPs. Moreover, a new scenario 
for VHE blazars emerges, wherein all values of r em (-z) fall within a small strip about the 
horizontal best-fit regression line in the r em — z plane, and the large scatter of the observed 
values of the slope arises from the large spread of the blazar redshifts. Remarkably, this 
issue can be settled in the near future not only by means of astrophysical VHE data from 
the upcoming CTA (Cherenkov Telescope Array), HAWC (High-Altitude Water Cherenkov 
Observatory), GAMMA-400 (Gamma Astronomical Multifunctional Modular Apparatus), 
LHAASO (Large High Altitude Air Shower Observatory) and HiSCORE (Hundred Square 
km Cosmic Origin Explorer), but also thanks to laboratory data based on the planned exper¬ 
iments ALPS II (Any Light Particle Search) and IAXO (International Axion Observatory). 
Our Universe may in this way be offering us a compelling reason to push physics beyond the 
Standard Model along a very specific direction, which possibly sheds light also on the nature 
of cold dark matter. 
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I. INTRODUCTION 

Nearly 10% of Active Galactic Nuclei (AGN) possess two oppositely-oriented collimated beams 
which emit photons up to the TeV energy range. When one of the beams is oriented towards us, 
the AGN is called a blazar. 

Thanks to the observations carried out with the Imaging Atmospheric Cherenkov Tele¬ 
scopes (IACTs) like H.E.S.S. [1], MAGIC [2] and VERITAS j3], according to the Tevcat cat¬ 
alog 43 blazars with known redshift have been detected in the very-high-energy (VHE) range 
(100 GeV- 100 TeV) gj. We stress that 40 of them are in the flaring state, whose typical lifetime 
ranges from a few hours to a few days. As it will explained below, 3 of them 1ES 0229+200, PKS 
1441+25 and S3 0218+35 will be discarded for the sake of the present analysis. 

All observed spectra of the considered VHE blazars are well fitted by a single power-law, and 
so they have the form d+bs^o, z) oc Kobs,o{z) E 0 ° , where Eq is the observed energy, while 

Kobs,o( z ) and T 0 b s (^) denote the normalization constant and the observed slope, respectively, for 
a source at redshift 2 . 

Unfortunately, the observational results do not provide any direct information about the intrinsic 
properties of the sources, as the VHE gamma-ray data strongly depend on the nature of photon 
propagation. Indeed, according to conventional physics the blazar spectra in the VHE range 
are strongly affected by the presence of the Extragalactic Background Light (EBL), namely the 
infrared/optical/ultraviolet background photons emitted by all galaxies since their birth (for a 
review, see 0)- However, it should be kept in mind that if some yet-to-be-discovered new physics 
changes the photon propagation, then some intrinsic source properties that are currently believed 
to be true may actually be incorrect. 

We restrict our discussion to the two standard competing models for VHE photon emission 
by blazars, namely the Synchrotron-Self-Compton (SSC) mechanism |0, [7j and the Hadronic Pion 
Production (HPP) in proton-proton scattering [8j. Both predict emitted spectra which, to a good 
approximation, have a single power-law behavior <h em (-E) = K em E~ r ° m for all the considered VHE 
blazars, where I \ em is the normalization constant and r em is the emitted slope. 

The relation between 4> 0 b s (-Eo, z) and & em (E) can be expressed in general terms as 

$ob S {E 0 ,z) = P 7 ^ 7 (E 0 ,+) $em(-Eo(l + z)) , (1) 

where P 7 ^. 7 (£'o, z) is the photon survival probability from the source to us, and is usually written 
in terms of the optical depth Tj(Eo,z) as 

P 1 ^(E 0 ,z)=e- T ^ Eo ’ z K ( 2 ) 

Before proceeding, a remark is compelling. A few years ago, a radically different mechanism 
has been put forward in order to explain the I ACT observations. Basically, the idea is that protons 
are accelerated inside blazars up to energies of order 10 11 GeV, while VHE emitted photons are 
neglected altogether. When the proton distance from the Galaxy is in the range 10 — 100 Mpc, 
they scatter off EBL photons through the process p + 7 —> p + 7 r°, so that the immediate decays 
vr° —>• 7 + 7 produce an electromagnetic shower of secondary photons: it is these photons that 
replace the emitted photons in such a scenario |]5]. Two characteristic features of this mechanism 
should be stressed in connection with the present analysis. 

• It is expected this effect to be important for sources at redshifts z > 0.15 and energies 
Eq > 1 TeV [10]. 

• Observed blazars variability shorter than 0.1 yr can be explained only for z > 0.20 at Eq > 
1 TeV HO]. 
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Now, a glance at the values of z for the considered sources reported in Table III shows that we have 
28 blazars with z < 0.15, and so - owing to the first item - these blazars cannot be explained by 
the mechanism in question. Let us next turn to the issue of variability. The overwhelming majority 
of the sources considered here are flaring, with a lifetime typically ranging from a few hours to a 
few days. Only the blazars 1ES 0229+200, 1ES 0806+524 and 1ES 1101-232 have a constant 7-ray 
luminosity ED- A look at the energy range AEq over which each source is observed (reported in 
Table III) shows that the lowest values of Eq are much below 1 TeV (with one exception to be 
addressed below). Because the shapes of the observed spectra do not exhibit any peculiar feature, 
it is evident that a single blazar is explained by the a single mechanism , and so we must conclude 
that the considered approach can explain at most 3 sources. Thus, here we do not address such an 
alternative possibility - whose relevance can become important at much larger energies - up to an 
exception. An analysis of the properties of the blazar 1ES 0229+200 has shown that it can hardly 
be explained by the SSC mechanism. Moreover, since its VHE luminosity is constant, this source 
is more likely to be explained by the proton emission model m- For this reason, we discard it 
from our discussion. 

The main issue we are concerned with in this paper is a possible correlation between the distri¬ 
bution of blazar VHE emitted spectra and the redshift. 

Superficially, the reader might well wonder about such a question. Why should a correlation of 
this kind exist? Cosmological evolutionary effects are certainly harmless out to redshift z — 0.5, 
and when observational selection biases are properly taken into account no such a correlation is 
expected to show up. As we shell see, this is not the case, and it is indeed this fact that has 
prompted our analysis. 

We are now in position to outline the structure of the paper. In Section II we report all 
observational information needed for the present analysis. Section III is devoted to inferring for each 
considered source the emitted slope r^(-z) and the emitted flux normalization (z) starting 
from the observed ones T 0 b s (z) and /t 0 bs,o(-2 ; ) assuming conventional physics, namely taking into 
account the effect of the EBL absorption alone. We next plot the values of T)+ ( z ) and ( z ) 

as a function of redshift z for all considered blazars. After performing a statistical analysis of 
the {r^'(z)} distribution, we end up with the conclusion that the resulting best-fit regression 
line decreases with increasing redshift. While this trend might be interpreted as an observational 
selection effect, a deeper scrutiny based on observational information shows that this is by no means 
the case. So, we are led to the conclusion that such a behavior is at odd with physical intuition, 
which would instead demand the best-fit regression line to be redshift-independent. In Section IV 
we try to fit the same data in the T^ — 2 plane with a horizontal straight line, but we conclude 
that it does not work. As an attempt to achieve a physically satisfactory scenario, in Section 
V we introduce Axion-Like Particles (ALP) (for a review, see (TSQ). They are spin-zero, neutral 
and extremely light pseudo-scalar bosons predicted by several extensions of the Standard Model 
of particle physics, and especially by those based on superstring theories [T3H23] . They interact 
only with two photons. Depending on their mass and two-photon coupling, they can be quite good 
candidates for cold dark matter [23] and give rise to very interesting astrophysical effects (to be 
discussed in Section VI), so that nowadays ALPs are attracting growing interest. Specifically, we 
suppose that photon-ALP oscillations take place in extragalactic magnetic fields of strength about 
0.1 nG - in agreement with the predictions of the galactic outflows models [251125] - as first proposed 
in m- As a consequence, photon propagation gets affected by EBL as well as by photon-ALP 
oscillations, whose combined effect is to substantially reduce the cosmic opacity brought about by 
the EBL alone, thereby considerably widening the conventional 7-ray horizon |28| . Accordingly, we 
re-derive for every source the emitted slope r^ p (z) and the emitted flux normalization K^ p (z) 
starting from the observed ones Tobs^) and K o ^ s0 (z). Proceeding as before, we plot again the 


values of r^ 1 pp (2:) and K^ F (z) as a function of redshift 2 for all considered sources. A statistical 


-ALP 
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analysis of the {r^^z)} distribution now shows that for a realistic choice of the parameters the 
corresponding best-fit regression line turns out to be indeed independent of redshift. Moreover, the 
values of r^ p (^) for the individual blazars exhibit a small scatter about such a best-fit straight 
regression line. We stress that the fact that the best-fit regression becomes a straight line horizontal 
in the — z plane is an amazing fact, which is the only one in agreement with physical intuition 
and provides a strong hint of the existence of ALPs. Finally, in Section VI we briefly discuss a 
new view of VHE blazars emerging from our result, its relevance - also in connection with other 
VHE astrophysical achievements employing ALPs - and its implications for the future of VHE 
astrophysics. Finally, in order to avoid breaking the main line of thought by somewhat involved 
technicalities concerning the evaluation of the photon survival probability in the presence of ALPs, 
we report this matter in the Appendix. 


II. OBSERVATIONAL INFORMATION 

The observational quantities concerning every blazar which are relevant for the present analysis 
are: the redshift z, the observed flux <h 0 b s (Eo, z) and the energy range AEo where each source is 
observed. 

Some care should be payed to the normalization constant entering the expression of <I> 0 bs(Eo, z), 
which is usually written as 


< f ) obs(A'o, ^0 — -^Q>bs,0C^) 


E 0 \ r ° ba(2) 

E re f J 


( 3 ) 


where E re { is an arbitrary energy value needed to make the exponential dimensionless, and in general 
varies from source to source. Manifestly, for the sake of comparison among the flux normalization 
constants for the different sources we have to recast Eq. ([3]) into a form such that E re f gets replaced 
by a quantity which is equal for all sources. Choosing as a fiducial normalization energy 300 GeV, 
it amounts to defining the new normalization constant as 


-^obsC^) — AT 0 bs,o(^) 


300 GeV\ _r ° ba(2) 


E, 


ref 


( 4 ) 


so 


that Eq. © becomes 


4 > obs(^'0)-2-) ^-obs ( z ) 


Ei 


0 


300 GeV 


TobsO) 


( 5 ) 


Note that we have <l> o bs(300 GeV, z) = K 0 ^ s (z). Henceforth, we shall deal exclusively with the 
values of K 0 \ )S (z). Therefore, we need to know both iL 0 b s (^) and r o b s (2) for any source. 


The values of z, r o b s (V), AEo and K 0 ^ s (z) for all considered blazars are reported in Table III 
Observe that while the error bars associated with r o b s (^) are quoted, those referring to K a \ )S (z) 
are not because of their lack in many of the published papers, but for our consideration this is not 
a big problem (nevertheless, the reader should keep this point in mind throughout the paper). The 
observed slope r o b s and the normalization constant K 0 bs as a function of z are plotted in the left 
and right panels of Fig. [lj respectively. 

Finally, denoting by V 0 b s the total number of detected photons per unit time, since <f> 0 bs(Eo, z) = 
dN^/(dEo dA) it follows that observed 7-ray luminosity per unit area E 0 b s , ae 0 ( z ) is the integral 
of $ ohs (E 0 ,z) over AE o- 
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FIG. 1: Left panel : The values of the slope T 0 bs are plotted versus the source redshift z for all considered 
blazars. Right panel: The values of the normalization constant A' 0 b s are similarly plotted versus the source 
redshift 0 for the same blazars in the left panel. 


III. CONVENTIONAL PROPAGATION IN EXTRAGALACTIC SPACE 


After a long period of uncertainty on the EBL precise spectral energy distribution and photon 
number density, today a convergence seems to be reached [5], well represented e.g. by the model 
of Franceschini, Rodighiero and Vaccari (FRV) [29], which we use for convenience. 

Owing to 77 —> e + e~ scattering off EBL photons mm, the emitted VHE photons undergo 
an energy-dependent absorption, so that the VHE photon survival probability is given by Eq. ([ 2 ]) 
with T<y(Eo,z) — > t rrv (Eo,z), where t rrv (Eo,z) is the optical depth of the EBL as evaluated 
within the FRV model in a standard fashion using the photon spectral number density |32h34| . As 
a consequence, the observed flux <f> 0 b s (Eo, z) is related to the emitted one (E) by 

d> obs (£b, z) = e-< HV ( E °’*) <h e c £(£ 0 (1 + z)) . (6) 


Let us begin by deriving the emitted spectrum of every source, starting from the observed one, 
within conventional physics. As a preliminary step, we rewrite Eq. © as 

( I 1 \ r obs (-^) 

30o4v) • (7) 

Owing to the presence of the exponential in the r.h.s. of Eq. ©, (Eq(1 + z)) cannot behave 

as an exact power law (unless t^ rv (Eo,z) has a pure logarithmic z-dependence). Yet, we have 
pointed out that it is expected to be close to it. So, we best-fit (£b(l + z)) hr Eq. © to the 

single power-law expression 


« BF (£o(i + *)) 



/ (1 + z)Eq 
V 300 GeV 


-rS« 


( 8 ) 


over the energy range AEo where the source is observed. Incidentally, we neglect error bars in 
the values of t rrv (Eo, 2 :) because they are not quoted by FRV. Since T^ R (z) depends linearly on 
r o bs(")> the inferred values of r^(z) have the same error bars as the values of r obs (z) quoted in 
Table [LEU 
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Of course, we are well aware that the best procedure would be to de-absorb each bin of the 
observed spectrum, thereby getting the emitted spectrum and next applying to it the above best¬ 
fitting procedure. Unfortunately, such a strategy is not viable in practice because the single ob¬ 
served energy bins with related error bars are not available from published papers. Nevertheless, 
the difference between the two procedures is expected to be relevant only for those sources whose 
highest energy points are affected by a large uncertainty, like 1ES 0229+200 - observed up to about 
11 TeV - which is however discarded. 

The values of the emitted slope T^ (z) are reported in Table IV and plotted in the le ft p anel 
of Fig. [ 2 J Similarly, the values of the normalization constant K ^ ( z ) are listed in Table VI and 
plotted in the right panel of Fig. [2j Again, denoting by N em the total number of emitted photons 
per unit time, because $ e m (E,z) = dN em /(dEdA), we see - similarly as before - that the emitted 
7 -ray luminosity per unit area F^ AE (z) is the integral of <f> em (F, z) over A E, whose values are 


reported in Table VII 
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FIG. 2: Left panel : The values of the slope T^J’ are plotted versus the source redshift 2 for all considered 
blazars. Right panel: The values of the normalization constant are similarly plotted versus the source 
redshift z for the same blazars in the left panel. 

We proceed by performing a statistical analysis of all values of T~j (z) as a function of z. Specifi¬ 
cally, we use the least square method and try to fit the data with one parameter (horizontal straight 
line), two parameters (first-order polynomial), and three parameters (second-order polynomial). In 
order to test the statistical significance of the fits we evaluate the corresponding Xred- The va l ues 
of the Xr e d obtained for the three fits are \^ ed = 2.35, y^ ed = 1.83 and x(? ed = 1.87, respectively. 
Thus, data appear to be best-fitted by the first-order polynomial 

T°£(z) = 2.68-2.21 2 . (9) 

The of {Tg^} distribution as a function of z and the associated best-fit straight regression line 
as defined by the last equation are plotted in Fig. [3j 

We stress that in order to appreciate the physical meaning of the best-fit straight regression line 
in question we should recall that (z) is the exponent of the emitted energy entering <h~^(.E). 
Hence, in the two extreme cases we have 


$eZ(E,0) oc E~ 2 - 68 


$^(E,0.6)ocE- h35 


( 10 ) 
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FIG. 3: Same as the left panel of Fig. [2] but with superimposed the best-fit straight regression line given 
by Eq. [9} 


thereby implying that its nonvanishing slope gives rise to a large variation of the emitted flux with 
redshift. 

Actually, one of the effects of the obtained best-fit straight regression line is that blazars with 
harder spectra are found only at larger redshift. What is the physical meaning of this fact? 

Since we have intentionally neglected the two blazars PKS 1441+25 and S3 0218+35 both at 
z — 0.94, our set of sources extends up to z ~ 0.54 (3C 279). Therefore, we are concerned with a 
relatively local sample, and so cosmological evolutionary effects are insignificant. 

Let us proceed to address all possible observational selection biases. 

• As we look at larger distances only the brighter sources are observed while the fainter ones 
progressively disappear. 

• Looking at greater distances entails that larger regions of space are probed, and so - under 
the assumption of an uniform source distribution - a larger number of brighter blazars should 
be detected. Now, the physical explanation of the best-fit straight regression line in Fig. [3] 
would naturally arise provided that F^ AE (z) tightly correlates with r^’(z) in such a way 
that brighter sources have harder spectra. Then the above selection bias translates into the 
statement that looking at greater distances implies that a larger number of blazars with 
harder spectra should be observed, which is just what Fig. [3] tells us. So, the real question 
concerns the existence of a tight correlation between AK (z) and rj+j’ (z). 

• Similarly, a logically consistent possibility would be that the jet opening angle 8(z) were 
tightly correlated with r ~^* (z) so that sources with stronger beaming had harder spectra: in 
such a situation - under the previous assumption of an uniform source distribution - the 
probability that the beam points towards us increases as larger regions of space are probed, 
and so sources with harder spectra would be more copiously found at larger distances, again 
in agreement with what Fig. [3] entails. 

In order to get a deeper insight into the issue addressed in the second item, we proceed by 
plotting F^ AE versus r£^* in Fig. |4j We see that F^ AE and F^’ are totally uncorrelated. As 
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FIG. 4: F^ ae is plotted versus 


a consequence, the cheap explanation outlined in the second item concerning the behavior of the 
best-fit straight regression line in Fig. [3] is doomed to failure. 

Finally, the possibility contemplated in the third item can also be excluded very easily by 
comparing the values of F^ (z) as reported in Table IV with those estimated for 8(z) in [55]. 


What turns out is that no correlation whatsoever between 8(z) and F^ (z) exists. 

It is very difficult to imagine an intrinsic mechanism which could explain the best-fit straight 
regression line in Fig. [3] in a physically satisfactory way within conventional physics. Otherwise 
stated, how can a source get to know its redshift z in such a way to adjust its emitted slope F^ (z) 
so as to reproduce the distribution with best-fit straight regression line reported in Fig. [3]? 

We are therefore led to regard the situation emerging from the above discussion as in manifest 
disagreement with physical intuition, which would instead demand the best-fit straight regression 
line to be horizontal in the — z plane. 


IV. AN ATTEMPT BASED ON CONVENTIONAL PHYSICS 

In spite of our previous finding, let us nevertheless try to impose by hand that the same data 
set {r^ > (z)} considered before is fitted by a horizontal straight line, and see what happens. The 
result is exhibited in Fig. [HJ In this case, we have = 2.41 and %[? ed = 2.35. 

Manifestly, this scenario does not work, since the value of \^ cd is unduly large. This fact hardly 
comes as a surprise, since we are not best-fitting the data. Still, this is a useful exercise, since 
it quantifies the price we have to pay in order to have a horizontal fitting straight line within 
conventional physics, and will be a benchmark for comparison when the ALP scenario will be 
considered. 

Much for the same reason, it is instructive to encompass 95 % of the observed sources inside a 
strip centered on the horizontal fitting line = 2.41. What is its width AT^? The answer is 
= 0.94, which is 39% of the value F^’ = 2.41. 
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FIG. 5: Horizontal fitting straight line in conventional physics. The values of r£j£ are plotted versus the 
source redshift z for all considered blazars with the corresponding error bars. Superimposed is the horizontal 
fitting straight line T^ = 2.41 with y^ ed = 2.35. The grey strip encompasses 95% of the sources and its 
width is = 0.94, which equals 39% of the value T^ = 2.41. 


V. AN ATTEMPT BASED ON AXION-LIKE PARTICLES 

As an alternative possibility to achieve a physically satisfactory situation, we invoke new physics 
in the form of axion-like particles (ALPs). As discussed in the Appendix, their most characteristic 
feature is to couple only to two photons with a coupling constant g a 77 according to the Feynman 
diagram shown in Fig. [ 6 j Clearly, in the presence of the extragalactic magnetic field B one photon 
line in Fig. [ 6 ] represents the B field, and so we see that in such a situation energy-conserving 
oscillations between VHE photons and ALPs take place m- Accordingly, photons acquire a split 
personality, traveling for some time as real photons - which suffer EBL absorption - and for some 
time as ALPs, which are unaffected by the EBL (as explicitly shown in the Appendix). As a 
consequence, t 7 (Eo,z) gets replaced by the effective optical depth t^ s (Eq,z), which is manifestly 
smaller than t 7 (Eo,z) and is a monotonically increasing function of Eq and z. The crux of the 
argument is that since the photon survival probability is now P^^(Eq, z) = e -T 7 ff l Eo ’ z ) ) even a 
small decrease of r 7 ff (Eo, z) with respect to rj( RV (i?o, z) gives rise to a large increase of the photon 
survival probability, as compared to the case of conventional physics. So, the main consequence 
of photon-ALP oscillations is to substantially attenuate the EBL absorption and consequently to 
considerably enlarging the conventional 7 -ray horizon [28] ■ 

Actually, P:^^ ( {Eq. z) can be computed exactly as outlined in the Appendix, where it is ex¬ 
plained that the extragalactic magnetic field likely has a domain-like structure (at least in first 
approximation), namely that B is homogeneous over a domain of size Lj om and has approximately 
the same strength B in all domains, but its direction randomly changes from one domain to the 
next [37H39J (more about this, in the Appendix). Furthermore, it will be shown that the considered 
ALP scenario contains only two free parameters £ oc g a77 B and Ld om : they obviously show up in 
Pj^y(Eo, z) even though such a dependence is not explicitly exhibited for notational simplicity. 
Realistic values of these parameters are £ = 0.1,0.5,1,5 and Ldom = 4Mpc, lOMpc, which will 
be regarded as our benchmark values. In reality, a third free parameter is the ALP mass m, but 
we assume m < 10~ 9 eV which entails that P^}^{Eq, z) is independent of m (as discussed in the 
Appendix). 





























11 



FIG. 6: Feynman diagram for the two-photon ALP coupling. 


At this point, knowing P^^(Eq, z), we proceed as above. To wit, we first write the emitted 


flux of every source as 


$ 


ALP 


(£0(1 + zj) = (p££(E 0 , zj) 1 K ohs (z ) 


En 


300 GeV 


Fobs (■ 2: ) 


(ID 


Next, we best-fit < h^' p (£'o(l + z)) in Eq. (11) to the single power-law expression 


^ALP,BF/ f /1 , „\\ _ F ALP /_\ f ( 1 + z ) E o\ rem ^ 

^ >em W + z))-K em (z) ^ 3Q0GeV j 


( 12 ) 


over the energy range AEo where the source is observed. It goes without saying that the remarks 
just below Eq. (J8j) apply here as well. 

This procedure is performed for each benchmark value of £ and Adom- We report in Table IY 
the values of T^ p (z) for Ad om = 4Mpc, £ = 0.1,0.5,1, 5, and similarly Table [V] contains the values 
of r^ pp (.s) for Tdom = lOMpc, £ = 0.1, 0.5,1, 5. 

We can at this point carry out a statistical analysis of the values of r^ p (2:) as a function of z, 
again for any benchmark value of £ and Adom- We still use the least square method and we try to fit 
the data with one parameter (horizontal line), two parameters (first-order polynomial) and three 
parameters (second-order polynomial). Finally, in order to quantify the statistical significance 
of each fit we compute the xf ed! whose values are reported in Table 0 for Adom = 4Mpc, £ = 
0.1, 0.5,1,5, and in Table [TT] for Ad om = 10 Mpc, £ = 0.1, 0.5,1, 5. 


# of fit parameters 

Xred.CP 

^red,ALP 



£ = 0.1 

£ = 0.5 

£ = 1 

£ = 5 

1 

2.35 

2.26 

1.43 

1.45 

1.55 

2 

1.83 

1.79 

1.46 

1.48 

1.57 

3 

1.87 

1.84 

1.48 

1.51 

1.60 


TABLE I: The values of Xred are displayed for the three fitting models considered in the text. In the 
first column the number of parameters in each fitting model is reported. The second column concerns 
deabsorption according to conventional physics, using the EBL model of FRV. The other columns pertain to 
the photon-ALP oscillation scenario with the EBL still described by the FRV model, and exhibit the values 
of xj? ed for Adom = 4 Mpc and different choices of £. The value in bold-face is the minimum of Xj? ed - 
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of fit parameters 

Xred, CP 

^red,ALP 



<T(\ 

II 

O 

l— 1 

£ = 0.5 

£ = 1 

£ = 5 

1 

2.35 

2.05 

1.39 

1.51 

1.55 

2 

1.83 

1.72 

1.43 

1.54 

1.57 

3 

1.87 

1.77 

1.47 

1.57 

1.60 


TABLE II: Same as Table [Ij but with L dom = lOMpc. The value in bold-face is the minimum of Xred- 


Hence, we see that the best-fitting procedure singles out the two following preferred cases. 


• L dom = 4Mpc, £ = 0.5, the horizontal straight regression line with equation r^ p = 2.52 
and Xred, ALP = L43 ‘ 

• L dom = lOMpc, £ = 0.5, the horizontal straight regression line with equation r^ p = 2.58 
and Xred,ALP = 


For simplicity, only for the two preferred cases are the values of r^ p (z) plotted in the left 
panels of Fig. [t} Also, only for these two cases are the values of K^ p (z) reported in Table VI and 
plotted in the right panels of Fig. [7j 

Clearly, for any choice of L& om the best-fitting procedure fixes the values of £ and y^ ed ALP . 
Even though we have taken realistic values for Ldom> in the absence of any information about its 
actual value it can well happen that other interesting results can emerge for different values of 
hdom' 


Let us summarize somewhat schematically the main achievements offered by the present sce¬ 
nario, especially emphasizing analogies and differences with respect to the situation based on 
conventional physics (considered in Sections III and IV). 


• The existence of an ALP with mass m < 10 -9 eV and suitable realistic values of the pa¬ 
rameters - such as £ = 0.5 and Ldom = 4Mpc or Ldom = lOMpc - gives rise to sizable 
photon-ALP oscillations in extragalactic space and yields a best-fit straight regression line 
for the {r^ P (£)} distribution which is just horizontal , in perfect agreement with the expec¬ 
tation based on physical intuition. We stress that this is an automatic consequence of the 
present model, and not an ad hoc requirement as in the case discussed in Section IV. The 
situation is illustrated in the left panels of Fig. [7j 


• The lack of correlation between F eni ^E and r em found in the context of conventional physics 
remains true in the presence of photon-ALP oscillations, thereby implying that again obser¬ 
vational selection biases play no role. This is shown in Fig. [9j 


VII shows that the values of are systematically slightly larger than 


A look at Table 

those of E c ^ p aa . This fact is in line with our expectations, and the reason is as follows. 
Since the values of -F 0 b s ,A.E are obviously fixed, such a difference has to be due to the photon 
propagation in the presence of photon-ALP oscillations and is the result of two competing 
effects. One of them is that at energies slightly larger than 100 GeV the EBL absorption 
is negligible. Accordingly, in the presence of many magnetic domains (as in the consid¬ 
ered situation) it is well known that the emitted photon flux gets reduced by a factor 1/3, 
due to the equipartition among the three degrees of freedom (2 photon polarization states 
and 1 ALP state), thereby producing a dimming of the source as compared to conventional 
physics [36]. The other effect occurs at larger energies, where EBL absorption becomes 
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important. Correspondingly, photon-ALP oscillation enhance the observed photon flux as 
compared to conventional physics, but since the emitted photon flux is larger al lower ener¬ 
gies, F^ J AE has to be larger than F^ AE in order to produce the same -Fobs.AB- 

• As in Section IV, it is interesting to compute the width Ar^ p of the strip - centered on the 
best-fit straight regression line r^ p - which encompasses 95% of the considered blazars. 
We find Ar^ p = 0.70 which is 28% of the value r^ pp = 2.52 for Ldom = 4Mpc, and 
AT^ PP = 0.60 which is 23% of the value r^ p = 2.58 for Ldom = lOMpc. These widths are 
considerably smaller than what we found within conventional physics in Section IV, namely 
Arf p = 0.94. 



FIG. 7: Left panels: The values of the slope T^ p are plotted versus the source redshift £ for all considered 
blazars. Right panels: The values of the normalization constant -fQm p are similarly plotted versus the source 
redshift 2 for the same blazars in the left panels. The upper row refers to the case £ = 0.5 and Tdom = 4Mpc, 
while the lower row to the case £ = 0.5 and Tdom = lOMpc. 
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FIG. 8: Same as the left panels of Fig. [7] but with superimposed the best-fit horizontal straight regression 
line. Moreover, in either case the grey band encompasses 95 % of the considered sources. 
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FIG. 9: Left panel: is plotted versus r^ p for £ = 0.5 and Ldom = 4Mpc. Right panel: Same as 

left panel but for £ = 0.5 and Ldom = lOMpc. 


VI. DISCUSSION AND CONCLUSIONS 

The main goal of the investigation reported in this paper concerns a possible correlation between 
the distribution of VHE blazar emitted spectra and the redshift. 

Broadly speaking, two logically distinct results have been obtained, which can somewhat 
schematically be summarize as follows. 

Working within conventional physics - and in particular within the two standard VHE photon 
emission models - we have shown that the emitted slope distribution {T^(2)} exhibits a correlation 
with z, since the associated best-fit regression line is a decreasing function of z given by T £^ (z) = 
2.68 — 2.21 z, with Xred = 1-83. This fact runs against physical intuition, since it is hard to 
understand how the source distribution can get to know the redshifts in such a way to adjust 
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their r em (z) values in such a way to reproduce such a statistical correlation. Indeed, this situation 
cannot be explained either by cosmological evolutionary effects or by observational selection biases. 
A further consequence is that blazars with harder spectra to be found only at larger redshift. 

We have shown that a way out of this conundrum involves new physics in the form of ALPs, 
with photon-ALP oscillations taking place in the extragalactic magnetic held. Their result is 
ultimately to substantially reduce the level of EBL absorption. We have focusses our attention 
on two realistic benchmark cases: Ldom = 4Mpc and Ldom = lOMpc. After having worked out 
the effect of photon-ALP oscillations on the VHE emitted spectra starting from the observed ones, 
we have discovered that in either case the best-fit regression line pertaining to the {r^ P (z)} 
distribution turns out to be a straight line in the r^ pp — z plane, hence independent of 2. This 
circumstance looks astonishing. Of course, by changing the effective level of EBL absorption we 
obviously expect the z-dependence of the best-fit {rsJf p (z)} distribution to differ from that of the 
( r em(A)} distribution. As a consequence, the inclination of the best-fit straight regression line in 
the T em — z plane should change. But to become exactly horizontal - which is the only possibility 
in agreement with physical expectation - looks almost like a miracle. In addition, the values of 
r^ pp (z) for the individual sources turn out to have a fairly small scatter about the considered 
best-fit straight regression line. 

Actually, the considered ALP framework possesses further nice features. 

• Even from a purely statistical point of view it is better than the conventional scenario, in 
which we have found x(? ed = 1.83 for the best-fit straight regression line (Section III) or 
A?ed = 2-35 for the horizontal fitting straight line (Section IV). Instead now the best-fit 
straight regression lines in the right panels of Fig. [7] have x'red = 1-43 for Ldom = 4Mpc and 
A red = 1-39 for Ldom = lOMpc. Note that this result has been obtained with a single free 
parameter. 

• A new picture arises wherein a sharp distinction between fundamental physics and boundary 
conditions naturally emerges. The above discussion implies that 95 % of the considered 
sources have a small spread in the values of r^ p (z). Specifically, r^ p (z) departs from the 
value of the best-fit straight regression line by at most 14% for £ = 0.5 and Ld om = 4Mpc 
and by 12% for £ = 0.5 and Ld om = lOMpc. Actually, the small scatter in the values of 
r^ p (z) implies that the emission mechanism is basically identical for all sources. On the 
other hand, the larger scatter in the values of K em (z) - presumably unaffected by photon- 
ALP oscillations when error bars are taken into account - is naturally traced back to the 
different environmental state of each source, such as for instance the accretion rate. A 
natural question should finally be addressed. How is it possible that the large scatter in 
{IWa)} distribution exhibited in the left panel of Fig. [I] arises from the small scatter in 
the {r^) PP (^)} distribution shown in the left panels of Fig. [%]? The answer is very simple: 
most of the scatter in the {Tobs^)} distribution arises from the large scatter in the source 
redshift. 

Before closing this Section, we find it worthwhile to put the result of our investigation into its 
proper perspective. 

During the last decade the interest in ALPs has been steadily growing. Various reasons have 
conspired towards this circumstance. With different motivations, the astrophysical implications 
of ALPs have been addressed since twenty years mmm- Certainly the claimed discovery of 
an ALP by the PVLAS collaboration [48] in 2005, which - even if subsequently withdrown by the 
PVLAS collaboration itself (45] - has provided a stimulus to look for astrophysical cross-checks ISB- 
[54]. Soon thereafter, it has been realized that several superstring theories predict the existence of 
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not only a single kind of ALP but more generally of a family of ALPs with very small masses nn- 
[25] , As a matter of fact, the study of the relevance of ALPs for high-energy astrophysics has 
continued until the present mmm- 

A tension between the predicted EBL level causing photon absorption and observations in the 
VHE range has been claimed from time to time (see e.g. [93,194]). but then a better determination 
of the EBL properties has shown that no problem exists. As already stressed at the beginning of 
Section III, nowadays the situation is different, since different techniques basically lead to the same 
EBL model. Yet, it has recently been claimed that VHE observations require an EBL level even 
lower than that predicted by the minimal EBL model normalized to the galaxy counts only [95], 
This is the so-called pair-production anomaly , which is based on the Kolmogorov test and so does 
not rely upon the estimated errors. It has thoroughly been quantified by a global statistical analysis 
of a large sample of observed blazars, showing that measurements in the regime of large optical 
depth deviate by 4.2 a from measurements in the optically thin regime Systematic effects 

have been shown to be insufficient to account for such the pair-production anomaly, which looks 
therefore real. Actually, the discovery of new blazars at large redshift like the observation of PKS 
1424+240 have strengthened the case for the pair-production anomaly [83]. Quite recently, the 
existence of the pair-production anomaly has been questioned by using a new EBL model and a y 2 
test, in which errors play instead an important role [96]. Because the Kolmogorov test looks more 
robust in that it avoids taking errors into account, we tend to believe that the pair-production 
anomaly is indeed at the level of 4.2 a. 

Amazingly, the existence of photon-ALP oscillations with the same realistic choice of the model 
parameters provides an excellent explanation for three completely different phenomena occurring 
in the VHE band. 

• The pair-production anomaly is naturally explained in terms of photon-ALP oscillations in 
extragalactic magnetic fields [TH E9l [83]. This should hardly come as surprise, since — as 
already explained - the ultimate effect of photon-ALP oscillation is to substantially lower 
the effective EBL absorption level. 

• According to conventional physics, flat spectrum radio quasars (FSRQs) should not emit 
above 20 — 30 GeV. This is due to the fact that higher energy photons accelerated in the jet 
enter - at a distance of about 1 kpc from the centre - the so-called broad-line region (BLR), 
whose high density of ultraviolet photons gives rise to an optical depth r ~ 15 owing to the 
same 77 —> e + e _ absorption process considered above. Even in this context, photon-ALP 
oscillations substantially lower the photon absorption level inside the BLR while still staying 
inside the standard blazar models, thereby allowing VHE photons to escape the BLR and 
be emitted, in remarkable quantitative agreement with observations m ■ We stress that the 
detection of VHE photons from FSRQs still represents a serious challenge for conventional 
models. 

• Our findings described in the present paper show that the conventional scenario of photon 
propagation in extragalactic space is seriously challenged. Nevertheless, a physically satisfac¬ 
tory picture emerges by considering photon-ALP oscillations in the extragalactic magnetic 
field. 

Altogether, such a situation evidently strongly suggests a preliminary evidence for the existence 
of an ALP, whose parameters make it a good candidate for cold dark matter [24]. Moreover, it looks 
tantalizing that the issue can definitively be settled not only with the advent of the new gamma-ray 
detectors like the CTA (Cherenkov Telescope Array) [97], HAWC (High-Altitude Water Cherenkov 
Observatory) [98] , GAMMA-400 (Gamma Astronomical Multifunctional Modular Apparatus) [99], 
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LHAASO (Large High Altitude Air Shower Observatory) [ 100] and HiSCORE (Hundred Square k m 
Cosmic Origin Explorer) [ TOT] , but also thanks to laboratory data based on the planned experiments 
ALPS II (Any Light Particle Search) |102) 1103] and IAXO (International Axion Observatory) 10 i . 
which will have the capability to discover an ALP with the properties assumed in the present 
analysis. Also experiments exploiting the techniques discussed in |105dl07] are quite promising for 
the laboratory detection of an ALP of this kind. 
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APPENDIX 

Here we summarize the main properties of the ALP scenario considered in the text and the 
evaluation of the photon survival probability P^^(Eq,z). 

A. General properties of axion-like particles (ALPs) 

As already stated, ALPs are spin-zero, neutral and extremely light pseudo-scalar bosons similar 
to the axion, namely the pseudo-Goldstone boson associated with the global Peccei-Quinn sym¬ 
metry U (1)pq proposed as a natural solution to the strong CP problem }108l . The axion interacts 
with fermions, two gluons and two photons, and one of its characteristic feature is the existence of 
a strict linear relationship between its mass m and the two-photon coupling constant g ail 

m = 0.7 k (jja~ n 10 10 GeVj eV , (13) 

where k is a model-dependent constant of order 1 [TO5] . Nevertheless, ALPs differ from axions in 
two respects: (1) their mass m and two-photon coupling constant g ail are unrelated , and (2) ALPs 
are supposed to interact only with two photons m • Hence, they are described by the Lagrangian 

£alp = ^ d M a - * m 2 a 2 - * g ail F IJM F lw a = * d fl a d^a - * m 2 a 2 + g ail E B a , (14) 

where a and m are the ALP field and mass, respectively, and E and B are the electric and magnetic 
components of the held strength F(F^ u is its dual). 

We shall be concerned with the particular case of a monochromatic photon/ALP beam emitted 
by a blazar at redshift z and traveling through the ionized extragalactic space, where a magnetic 
held B is supposed to exists with the domain-like structure described in Section VI. Therefore E is 
the electric held of a beam photon while B is the extragalactic magnetic held. In such a situation, 
the mass matrix of the photon-ALP system is non-diagonal, so that the propagation eigenstates 
differ from the interaction eigenstates. This fact gives rise to photon-ALP oscillations in the 
beam mm, much in the same way as it happens in a beam of massive neutrinos of different 
flavor (however for photon-ALP oscillations an external magnetic held is needed to compensate for 
the spin mismatch). In order to avoid notational confusion, the symbol E will henceforth denote 
the energy of the beam particles. 

Since we suppose that E m (more about this, later), we can employ the short wavelength 
approximation. As a consequence, the photon/ALP beam propagation turns out to be described 
by the following hrst-order Schrodinger-like equation with the time replaced by the y-coordinate 
along the beam m 

(i-jt + E + m\ il>(y) = 0 (15) 

with wave function 

( A x (y) \ 

V’(y) = A z (y) , (16) 

V a(y) J 

where A x (y) and A z (y) denote the photon amplitudes with polarization along the x- and z-axis, 
respectively, while a(y) is the amplitude associated with the ALP. It is useful to introduce the 
3-dimensional basis vectors {|7a;), \lz), |a)} where |7 X ) = (1,0, 0) T and \^ z ) = (0,1,0)^ represent 
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the two photon linear polarization states along the x- and 2-axis, respectively, and |a) = (0,0, 1) T 
denotes the ALP state. Accordingly, we can write if(y) as 

V’(y) = A x (y) |7 X ) + A z (y) \^ z ) + a(y) |o) . (17) 

Note that the quantity 7~L = -(E-\-M) plays formally the role of the Hamiltonian. Denoting by 


U(y,y 0 ) the transfer matrix - namely the solution of Eq. (15) with initial condition U(yo,yo) = 1 


- the propagation of a generic wave function across the domain in question can be represented as 

V’(y) =u(y,yo)^(yo) • (18) 

So far our attention has been restricted to the wave function, which describes a polarized beam. 
However, since in the gamma-ray band the polarization cannot be measured, we are forced to 
employ the density matrix p(y). As in quantum mechanics, it obeys the analog of the Von Neumann 


equation associated with Eq. (15), whose form is 


i^f- = pM^ — Mp . 

ay 


(19) 


As we shall see, EBL absorption implies that T~L - and so M and p(y) - are not self-adjoint. Hence 
U(y,yo) fails to be unitary. Nevertheless, it is trivial to check that the familiar relation 


p(y) =u(y,yo) p{yo)u J (y,yo) 


( 20 ) 


retains its validity. Moreover, the probability that the beam in the state pi at y$ will be found in 
the state p 2 at y is still given by the standard relation 


P, 


P1-»P2 


(yo,y) = Tr(p 2 ^(y,yo)Myo)^ T (y,2/o)) 


( 21 ) 


where it is supposed as usual that Trpi = Trp2 = 1. 

We stress that the advantage arising from the short wavelength approximation is that the beam 
propagation can formally be described as a three-level non-relativistic decaying quantum system. 

It is quite enlightening to start by considering the beam propagation over a single magnetic 
domain and to neglect EBL absorption [57) 111 ] (this attitude indeed makes sense, since EBL 
absorption is independent of photon-ALP oscillations). Because then B is homogeneous, we can 


choose the 2-axis along B so that B x = 0. Correspondingly, the matrix A4 entering Eq. (15) takes 
the form 


M = 


-u‘J2E 0 0 

0 — Co'p 1 /2^7 g ajy B/2 

0 g ail B/2 —m 2 /2E 


( 22 ) 


where w p i is the plasma frequency of the ionized intergalactic medium (IGM), which is related 
to the mean electron density n e by w p i = 3.69 • 10 _11 (n e /cm _3 ) 1 ' /2 eV (111) . So, by defining the 
oscillation wave number as 


^osc — 


2 E 


v (5a77 b^ 


1/2 


(23) 


it is a simple exercise to show that the photon-ALP conversion probability across the considered 
domain is 


R 


7 —>a 


=(^f) : 


Sill 


^osc -^do 


(24) 
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Thus, defining the energy threshold m 


e l = 


\m 2 — Wpj| 

2 Qa'y'y B 


(25) 


we see that in the strong-mixing regime - namely when condition E > El is met - Eq. (23) reduces 
to 


^osc — Sa77 B 


(26) 


and consequently the photon-ALP conversion probability becomes maximal as well as energy- 
independent. Finally, since B and g a 77 enter £alp only in the combination g a77 B, it proves 
instrumental to define the quantity 


£— (hg) (Sa'y-yW GeV) 


in terms of which Eq. (25) can be rewritten as 

25.64 


E l ^ 


m \ 2 
IQ” 10 eV > 


w P i 

IQ" 10 eV 


GeV . 


(27) 


(28) 


B. Model parameters 

Before proceeding further, it is worthwhile to discuss the allowed values of the free parameters 
of the model. In this field, a signal is considered a discovery if its statistical significance is at 
least 5cr. Accordingly, the bounds should be divided into two classes, depending on whether their 
statistical significance is larger or smaller than 5a. 

• Those with statistical significance larger than 5a prevent any discovery of an ALP, and to 
date only one bound belongs to this class: the one obtained by the CAST experiment at 
CERN, which reads g a 77 < 8.8 • ICE 11 GeV for m < 0.02eV [ 112] . 

• All other weaker bounds provide useful information, but still allow for the discovery of an 
ALP. Various bounds of this sort have ben derived, but we discuss here only those that we 
believe to be most relevant. One is based on the study of 39 Galactic globular clusters and 
slightly improves the CAST constraint giving g a ~ n < 6.6 • 10~ n GeV at the 2a level (the 
mass range of its validity is not specified) [86]. Another has been obtained by the H.E.S.S. 
collaboration, looking at the source PKS 2155 — 304, from the absence of the characteristic 
fluctuating behavior of the realizations of the beam propagation around El (more about 
this, later) and reads g 0 - r/ < 2.1 • 10 -11 GeV for 1.5 • 10~ 8 eV < m < 6.0 • 10~ 8 eV at 
the 2a level [80]. A similar strategy has been applied to the source to the Hydra galaxy 
cluster, yielding g ail < 8.3 • 10 -12 GeV for eV < rn < 7.0 • 10 _12 eV at the 2a level [STJ. 
Finally, the explosion of the 1987A supernova has been used to set an upper bound on g ail 
depending on the value of m. The logic is as follows. Along with a neutrino burst lasting 
about 10 seconds, also an ALP burst should be emitted during the collapse. Inside the 
Galaxy magnetic field, some of the ALPs should convert into X-rays and next be detected 
by the Solar Maximal Mission (SMM) satellite. From the absence of detection, originally 
the bounds g ail < 1.0 • 10 -11 GeV for m < 1.0 • 10 -9 eV |4l] and g ail < 3.0 • 10~ 12 GeV for 
m < 1.0 • 10 -9 eV have been set [221 • A very recent analysis updating these two early works 
has been carried out, especially using the present state-of-the art knowledge of supernova 
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explosions (even if we still do not know why they actually explode). However, an unavoidable 
uncertainty affects the detection efficiency of the SMM, for instance the systematic error 
and the effective area, even because SN1987A was observed perpendicularly to the viewing 
direction of the SMM. The authors do not quote the statistical significance of this bound, 
which is certainly less than 2 a [92j. 


As far as the extragalactic magnetic held B is concerned, the situation is much less clear-cut. 
A general consensus exists that it possesses a domain-like morphology (at least in first approxima¬ 
tion) [3TH39] . As mentioned in Section IV, it is supposed that B is homogeneous over a domain of 
size Tdom equal to its coherence length, with B randomly changing its direction from one domain 
to another but keeping approximately the same strength. Unfortunately, the values of both B and 
Adom are largely unknown. We assume B to be pretty large according to the predictions of the 
galactic outflows models [25, 26|. Actually, within these models it looks natural to suppose that 
Adom is of the order of the galactic correlation length, which amounts to take Ad om = (1 — 10) Mpc 
in the present Universe (z = 0), even though larger values cannot be excluded. Then Fig. 4 of |113j 
implies that B < 1 nG, and so by taking the most restrictive bound on g a 77 it follows from Eq. (27) 
that £ < 8. The upper bound on the mean diffuse extragalactic electron density n e < 2.7T0 -7 cm -3 
is provided by the WMAP measurement of the baryon density m, which - thanks to the above 
relation - translates into the upper bound u v \ < 1.92 • 10~ 14 eV. What about ml In order to 
maximize the effect of photon-ALP oscillations we work throughout within the strong mixing 
regime. Recalling that the lower end of the VHE band is E = 100 GeV, we then have to require 
El < 100 GeV. Putting everything together and using Eq. (28) we find m < 1.97 • lO -10 ^ 1 / 2 , and 
at any rate we should have m < 5 • 10 _1 ° eV. However, we remark that the size of the effect remains 
basically unchanged if the bound El < 100 GeV becomes less restrictive by a factor of a few, and 
so for all practical purposes we can assume m < 10” 9 eV. Finally, by combining the CAST bound 
g ail < 8.8 • 10 -11 GeV with the interaction term in £alp> it is straightforward to get the following 
order-of-magnitude estimate for the cross-sections <7(07 —> / + /~) ~ a(af ± -» 7 f±) < 10~ 5 °cm 2 
(here / denotes any charged fermion), which show that effectively ALPs do not interact with 
anything and in particular with the EBL. 


C. Strategy 

One point should now be emphasized. Due to the nature of the extragalactic magnetic field, 
the angle of B in each domain with a fixed fiducial direction equal for all domains (which we 
identify with the z-axis) is a random variable, and so the propagation of the photon/ALP beam 
becomes a A^-dimensional stochastic process , where denotes the total number of magnetic 
domains crossed by the beam. Therefore we identify the photon survival probability with its 
average value. Moreover, we shall see that the whole photon/ALP beam propagation can be 
recovered by iterating A7 times the propagation over a single magnetic domain, changing each 
time the value of the random angle. Thanks to the fact that B is homogeneous in every domain, 
the beam propagation equation can be solved exactly in a single domain. Clearly, at the end we 
have to average the photon survival probability as evaluated for one arbitrary realization of the 
whole propagation process - namely for a particular choice of the random angle in each domain - 
over all possible realizations of the considered stochastic process (i.e. over all values of the random 
angle in each of the Nj domains) [ 115] . 

Our discussion is framed within the standard ACDM cosmological model with Dm = 0.3 and 
Ha = 0.7, and so the redshift is the natural parameter to express distances. In particular, the 
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proper length L^ ora {z a , Zb) extending over the redshift interval [z a ,Zb] is 


L(z a ,z h )~ 4.29-10 3 


<"Zb 


dz 


(1 + z)[0.7 + 0.3(1 + z) 3 ] 1 / 2 


Mpc ~ 2.96-10 3 In ( - — ) Mpc . (29) 

\l + lA5z a J F V ' 


Accordingly, the overall structure of the cellular configuration of the extragalactic magnetic field 
is naturally described by a uniform mesh in redshift space with elementary step A z, which is 
therefore the same for all domains. This mesh can be constructed as follows. We denote by 
Lj’jn = L({n — l)Az,nAz) the proper length along the y-direction of the generic n-th domain, 
with 1 < n < Nj, where the total number of magnetic domains Nj towards the considered blazar 
is the maximal integer contained in the number z/ Az, hence — z/Az. In order to fix Az we 
consider the domain closest to us, labelled by 1 and - with the help of Eq. (29) - we write its 


proper length as (L^ m /5Mpc) 5 Mpc = L(0, Az) = 2.96 • 10 3 ln(l + 1.45 Az) Mpc, from which 

we get Az ~ 1.17 • 10 -3 (L^ m /5 Mpc). So, once Lfo m is chosen in agreement with the previous 
considerations, the size of all magnetic domains in redshift space is fixed. At this point, two further 
quantities can be determined. First, the total number of the considered domains is Nj ~ zjAz — 
0.85 • 10 3 (5Mp C /4 1 o ) Jz. Second, the proper length of the n-th domain along the y-direction 
follows from Eq. (|29[) with z a — > (n — l)Az, Zf, —» n Az. Whence 


4oin^ 2.96-10 3 ln 


1 + 


1.45 Az 


1 + 1.45 (n- l)Az 


Mpc 


(30) 


D. Propagation over a single domain 

What still has to be done is to take EBL absorption into account and to determine the magnetic 
field strength in the generic n-th domain. 

The first goal can be achieved by simply following the discussion of [115] . Because the domain 
size is so small as compared to the cosmological standards, we can safely drop cosmological evo¬ 
lutionary effects when considering a single domain. Then as far as absorption is concerned what 
matters is the mean free path A 7 for the reaction 77 —> e + e _ , and the term i/ 2A 7 should be 
inserted into the 11 and 22 entries of the M matrix. In order to evaluate A 7 , we imagine that two 
hypothetical sources located at both edges of the n-th domain are observed. Therefore, we insert 
Eq. Q into Eq. (jTj) of Section I and further apply the resulting equation to both of them. With 
the notational simplifications < 3? 0 b s (-E'o, z) —> $(Eq) and 4> em (£'o(l + z)j -+ ^(^(l + z)), we have 

&(Eq) = e -+4o,(n-i)A 2 ) $(# 0 [i + ( n _i)A z\) , &(Eo) = e~ T ^ Eo ’ nAz) <5>(E 0 {l + nAz)) , (31) 

which upon combination imply that the flux change across the domain in question is 

<h(£ 0 [l + (n— l)Az]) = e -hW,«A*K(^,MA*)] $(£ 0 (l+nAz)) . (32) 

Because cosmological effects can be neglected across a single domain, Eq. ( f32| ) should have the 
usual form 

T(£o[ 1 + (n - l)Az]) = e _L do ) m/ A 7 l) (- E o) (1 + n Az)) , (33) 

and the comparison with Eq. ( |32| ) ultimately yields 

rW 

_^dom_ 

^(EqjU Az) — Tj(Eq, {n — l)Az) 


a4(E 0 ) 


(34) 
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where the optical depth is evaluated as stated in the text for the EBL model of FRV 

To accomplish the second task, we note that because of the high conductivity of the IGM medium 
the magnetic flux lines can be thought as frozen inside it [39]. Therefore the flux conservation during 
the cosmic expansion entails that B scales like (1 + z) 2 , so that the magnetic field strength in a 
domain at redshift 2 is B(z) = B(z = 0)(1 + z) 2 . Hence in the n-th magnetic domain we have 
B (n) = £(!)(! + ( re _i)Az) 2 . 


So, at this stage the matrix A4 in Eq. (15) as explicitly written in the n-th domain reads 


/ 


M (n) = 


i/ 2 A. 

0 


(n) 

‘7 


i/2 A. 


(n) 


\ sin V’n 5077/2 B ( n ) cos ip n g ail /2 


B (n) sin tp n g ail /2 \ 
B^ COSl/jngayry/2 

0 J 


(35) 


where ip n is the random angle between Bl n l and the fixed fiducial direction along the z-axis (note 
that indeed MP y A4). Apart from ip n , all other matrix elements entering are known. 

Finding the transfer matrix of Eq. (15) with given by Eq. (35) is a straightforward even if 

somewhat boring game. The result is 


Un(E n ,i/> n ) = e lEnL d°™ 


(n) 


? V 1 Ti(ip n ) + e v 2 d °W T 2 (ip n ) + e V 3 T 3 (ip n 


(36) 


with 


cos 2 ip n — sin ip n cos ip n 0 
T\(ipn) = | - sin ip n cos ip n sin 2 Ipn 0 

0 0 0 


(37) 


/ 


T 2 (lpn) = 


1 +V 1 • 2 /, 

8111 ^ 

-1+JI-A&1 . . . 

- / sin w n cos w r 

2^/1-452 Yn 


—14-v/l—4<52 , , ,•* . , 

- / sm w n cos w n 7 » sin w, 

-l+y / l-4<52 2 ./. ii5„ _/. 


\ 


- / COS 2 lp r 


iSn. 


\ \/i-4«S 


sin ipn 


7 =, g COS 1p n 
V l-45n 2 


," )n COS Ip. 
i+yi-4^ 

V 1 - 452 / 


(38) 


/ 1+a/1—452 . 2 , 1+a/ 1—452 . , , _j<5 . , \ 

- Y h Sill 2 Ipn - y A sm Ipn COS Ipn , Wn „ Sin Ipn 

" 2^/1^452 V ' 71 ^ 


r 3 (^„) = 


yi-< 


l+\/l —452 

V 1 - 452 


\ \/I Z 452 


sin 'Ipn COS Ipn 
sin ipn 




COS 


iSn 


1+a/ 1-452 ^ o 2 ,, 


COS Ipn 


~H"\/ l~ 4( ^n 
2\/l-45 3 / 


(39) 


where we have set 


A< n) = 


2A^ ) (£ 0 ) 


x M _ 

a 2 — 


4 A. 


(n) 


(1 - yn^ 2 ) , a^ } = (1 + yr^ 2 ) (40) 

4 A,, 


with 


E n — Eq 


1 + (n — 1) Az) 


d n ^U^\E 0 ) 


nG 


10 11 GeV 


( 41 ) 


where £ n is just £ as defined by Eq. (27) and evaluated in the n-th domain. 
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E. Calculation of the photon survival probability in the presence of photon-ALP oscillations 


Our aim is to derive the photon survival probability P^ pp {Eq, z) entering the text. So far, we 
have dealt with a single magnetic domain but now we enlarge our view so as to encompass the 
whole propagation process of the beam from the source to us. This goal is trivially achieved thanks 
to the analogy with non-relativistic quantum mechanics, according to which — for a fixed arbitrary 
choice of the angles <n<N d ~ the whole transfer matrix describing the propagation of the 

photon/ALP beam from the source at redshift z to us is 


N d 


U {Eq, z\ ip\, ...,ip Nd ) = U n (E n , ip n ) 


(42) 


71=1 


Moreover, the probability that a photon/ALP beam emitted by a blazar at z in the state p± will 
be detected in the state p 2 for the above choice of {i/j n }i<n<N d is given by Eq. plj). Whence 


p Pi^P2 ( E o,z-,ipi, -,i>N d ) = Tr (/ o 2 U {E 0 ,z-,il>i,...,i ’ Nd ) P\U X {E 0 ,z-,ipi,-A nS) ( 43 ) 
with Trpi = Trpa = 1- 

Since the actual values of the angles {i> n }i<n<N d are unknown, the best that we can do is to 


evaluate the probability entering Eq. (43) as averaged over all possible values of the considered 
angles, namely 


P p\ — yp2 (Eq, z') /Pp 1 —^p2 {Eq, z, Ipl, ..., V’iV’d) / 

\ / Pi,... 


p 1 ,...,pN d 


(44) 


indeed in accordance with the strategy outlined above. In practice, this is accomplished by evalu¬ 


ating the r.h.s. of Eq. (43) over a very large number of realizations of the propagation process (we 
take 5000 realizations) randomly choosing the values of all angles {^n}i<n<N d for every realization, 
adding the results and dividing by the number of realizations. 

Because the photon polarization cannot be measured at the considered energies, we have to 
sum the result over the two final polarization states 


1 0 0 
Px= | 0 0 0 
0 0 0 


(45) 


0 0 0 
Pz= I 0 1 0 
0 0 0 


(46) 


Moreover, we suppose here that the emitted beam consists 100 % of unpolarized photons, so that 
the initial beam state is described by the density matrix 


1 


1 0 0 


Punpol — q I 0 10 
' 0 0 0 


We find in this way the photon survival probability P^ p {Eq, z ) 


(47) 


P ALP {E 0 ,z) = /Pp unpol ^p x {Eq, z; ,..., i/> Nd ) ) + 

\ / pi ,...,pN d 

+ ( P Pun P ol^Pz {EOl Z\ Vh, 1> Nd ) ) 

\ / pl,...,p N 


( 48 ) 
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which is indeed the quantity entering Eq. © of the text. 

A final remark is in order. It is obvious that the beam follows a single realization of the 
considered stochastic process at once, but since we do not know which one is actually selected 
the best we can do is to evaluate the average photon survival probability. Nevertheless, it can be 
shown that the considered realizations exhibit an oscillatory behavior for an energy close to the 
threshold El defined by Eq. (25). Therefore, observations performed at energies close enough to 
El should exhibit oscillations in the energy spectrum [73j[82]. This fact has been used to derive 
the previous bounds discussed in f80 ; , 8 lj. 
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Source 

z 

r obs 

AE 0 [TeV] 

AVbs [cm 

~ 2 s _i TeV -1 ] 

3C 66B 

0.0215 

3.10 + 0.37 

0.12- 1.8 

1.74 

• io- 11 

Mrk 421 

0.031 

2.20 + 0.22 

0.13-2.7 

6.43 

• 10~ 10 

Mrk 501 

0.034 

2.72 + 0.18 

0.21 - 2.5 

1.53 

• IO" 10 

1ES 2344+514 

0.044 

2.95 + 0.23 

0.17-4.0 

5.42 

• 10- 11 

Mrk 180 

0.045 

3.30 + 0.73 

0.18-1.3 

4.50 

• 10- 11 

1ES 1959+650 

0.048 

2.72 + 0.24 

0.19-1.5 

8.99 

• 10- 11 

1ES 1959+650 

0.048 

2.58 + 0.27 

0.19-2.4 

6.03 

• 10- 11 

AP LIB 

0.049 

2.50 + 0.22 

0.30-3.0 


7 

1ES 1727+502 

0.055 

2.70 + 0.54 

0.10-0.6 

9.60 

•10- 12 

PKS 0548-322 

0.069 

2.86 + 0.35 

0.32-3.5 

1.10 

• 10- 11 

BL Lacertae 

0.069 

3.60 + 0.43 

0.15-0.7 

5.80 

• 10- 10 

PKS 2005-489 

0.071 

3.20 + 0.19 

0.32 - 3.3 

3.44 

• 10" 11 

RGB J0152+017 

0.08 

2.95 + 0.41 

0.32-3.0 

1.99 

• 10- 11 

1ES 1741+196 

0.083 

? 

? 


7 

SHBL J001355.9-185406 

0.095 

3.40 + 0.54 

0.42 - 2.0 

7.05 

•10" 12 

W Comae 

0.102 

3.81 ± 0.49 

0.27- 1.1 

5.98 

• 10" 11 

1ES 1312-423 

0.105 

2.85 + 0.51 

0.36-4.0 

5.85 

• 10- 12 

VER J0521+211 

0.108 

3.44 + 0.36 

0.22- 1.1 

5.35 

• 10- 11 

PKS 2155-304 

0.116 

3.53 + 0.12 

0.21-4.1 

1.27 

• 10- 10 

B3 2247+381 

0.1187 

3.20 + 0.71 

0.15-0.84 

1.40 

• 10- 11 

RGB J0710+591 

0.125 

2.69 + 0.33 

0.37-3.4 

1.49 

• 10- 11 

H 1426+428 

0.129 

3.55 + 0.49 

0.28 - 0.43 

1.46 

• 10- 10 

1ES 1215+303 

0.13 

3.60 + 0.50 

0.30 - 0.85 

2.30 

• 10- 11 

1ES 1215+303 

0.13 

2.96 + 0.21 

0.095 - 1.3 

2.27 

• 10- 11 

RX J1136.5+6737 

0.1342 

? 

? 


7 

1ES 0806+524 

0.138 

3.60+1.04 

0.32-0.63 

1.92 

•10- 11 

1ES 0229+200 

0.14 

2.50 + 0.21 

0.58-11 

1.12 

• 10- 11 

1RXS J101015.9-311909 

0.142639 

3.08 + 0.47 

0.26-2.2 

7.63 

• 10- 12 

H 2356-309 

0.165 

3.09 + 0.26 

0.22-0.9 

1.24 

• 10- 11 

RX J0648.7+1516 

0.179 

4.40 + 0.85 

0.21 - 0.47 

2.30 

• 10- 11 

1ES 1218+304 

0.182 

3.08 + 0.39 

0.18- 1.4 

3.62 

• 10- 11 

1ES 1101-232 

0.186 

2.94 + 0.22 

0.28-3.2 

1.94 

• 10- 11 

1ES 0347-121 

0.188 

3.10 + 0.25 

0.30-3.0 

1.89 

• 10“ n 

RBS 0413 

0.19 

3.18 + 0.74 

0.30 - 0.85 

1.38 

• 10" 11 

RBS 0723 

0.198 

? 

? 


? 

1ES 1011+496 

0.212 

4.00 + 0.54 

0.16-0.6 

3.95 

•10- 11 

MS 1221.8+2452 

0.218 

? 

? 


? 

PKS 0301-243 

0.2657 

4.60 + 0.73 

0.25-0.52 

8.56 

•IO” 12 

1ES 0414+009 

0.287 

3.45 + 0.32 

0.18-1.1 

6.03 

• 10~ 12 

S5 0716+714 

0.31 

3.45 + 0.58 

0.18-0.68 

1.40 

• 10- 10 

1ES 0502+675 

0.341 

3.92 + 0.36 

? 


7 

PKS 1510-089 

0.361 

5.40 + 0.76 

0.14-0.32 

6.97 

•IO” 12 

3C 66A 

0.41 

4.10 + 0.72 

0.23 - 0.47 

4.00 

• 10- 11 

PKS 1222+216 

0.432 

3.75 + 0.34 

0.08-0.36 

1.71 

• 10- 10 

1ES 0647+250 

0.45 

? 

? 


7 

PG 1553+113 

0.5 

4.50 + 0.32 

0.23- 1.1 

4.68 

•10- 11 

3C 279 

0.5362 

4.10 + 0.73 

0.08 - 0.46 

9.86 

• 10- 11 

PKS 1424+240 

> 0.6035 

3.80 + 0.58 

0.14-0.5 

1.09 

• 10~ n 


TABLE III: Considered VHE blazars with known energy, redshift, spectral slope T 0 b s , energy range and 
normalization constant AT 0 b s - Statistical and systematic errors are added in quadrature to produce the total 
error reported on the measured spectral slope. When only statistical errors are quoted, systematic errors 
are taken to be 0.1 for H.E.S.S., 0.15 for VERITAS, and 0.2 for MAGIC. Sources with question marks lack 
information to perform our analysis and are discarded. For PKS 1424+240 only a lower limit on the redshift 
exists in the literature: thus, it is neglected in our analysis. 
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Source 

fCF 

L em 

pALP 

1 em 



£ = o.i 

£ = 0.5 

£ = i 

£ = 5 

3C 66B 

3.00 

3.00 

3.00 

3.01 

3.03 

Mrk 421 

2.05 

2.05 

2.05 

2.06 

2.10 

Mrk 501 

2.54 

2.54 

2.54 

2.55 

2.59 

1ES 2344+514 

2.71 

2.71 

2.71 

2.73 

2.79 

Mrk 180 

3.07 

3.07 

3.07 

3.09 

3.14 

1ES 1959+650 

2.46 

2.46 

2.47 

2.49 

2.55 

1ES 1959+650 

2.32 

2.32 

2.32 

2.35 

2.40 

AP LIB 

2.21 

2.21 

2.22 

2.24 

2.31 

1ES 1727+502 

2.52 

2.52 

2.52 

2.54 

2.58 

PKS 0548-322 

2.44 

2.44 

2.45 

2.51 

2.58 

BL Lacertae 

3.29 

3.29 

3.30 

3.33 

3.39 

PKS 2005-489 

2.77 

2.77 

2.79 

2.85 

2.91 

RGB J0152+017 

2.47 

2.47 

2.48 

2.56 

2.63 

1ES 1741+196 

? 

? 

? 

? 

7 

SHBL J001355.9-185406 

2.81 

2.81 

2.84 

2.94 

3.01 

W Comae 

3.17 

3.17 

3.20 

3.31 

3.38 

1ES 1312-423 

2.17 

2.17 

2.21 

2.34 

2.40 

VER J0521+211 

2.79 

2.79 

2.82 

2.93 

3.00 

PKS 2155-304 

2.81 

2.81 

2.86 

3.00 

3.05 

B3 2247+381 

2.60 

2.60 

2.64 

2.74 

2.80 

RGB J0710+591 

1.89 

1.89 

1.96 

2.12 

2.16 

H 1426+428 

2.85 

2.85 

2.89 

3.01 

3.08 

1ES 1215+303 

2.75 

2.75 

2.81 

2.97 

3.03 

1ES 1215+303 

2.35 

2.35 

2.40 

2.51 

2.55 

RX J1136.5+6737 

? 

? 

? 

? 

7 

1ES 0806+524 

2.70 

2.70 

2.77 

2.93 

3.00 

1ES 0229+200 

0.61 

0.62 

1.20 

1.43 

1.26 

1RXS J101015.9-311909 

2.20 

2.20 

2.29 

2.45 

2.49 

H 2356-309 

2.05 

2.05 

2.17 

2.35 

2.40 

RX J0648.7+1516 

3.45 

3.45 

3.55 

3.71 

3.77 

1ES 1218+304 

1.97 

1.97 

2.13 

2.31 

2.34 

1ES 1101-232 

1.72 

1.73 

1.96 

2.13 

2.13 

1ES 0347-121 

1.87 

1.87 

2.11 

2.28 

2.28 

RBS 0413 

1.88 

1.89 

2.07 

2.28 

2.32 

RBS 0723 

? 

? 

? 

7 

7 

1ES 1011+496 

2.90 

2.90 

3.06 

3.22 

3.26 

MS 1221.8+2452 

? 

? 

? 

7 

7 

PKS 0301-243 

2.93 

2.93 

3.27 

3.46 

3.49 

1ES 0414+009 

1.65 

1.65 

2.12 

2.26 

2.25 

S5 0716+714 

1.60 

1.61 

2.07 

2.22 

2.22 

1ES 0502+675 

? 

? 

? 

7 

7 

PKS 1510-089 

4.02 

4.02 

4.33 

4.45 

4.48 

3C 66A 

1.53 

1.54 

2.31 

2.40 

2.39 

PKS 1222+216 

2.46 

2.46 

2.79 

2.87 

2.89 

1ES 0647+250 

? 

? 

? 

7 

7 

PG 1553+113 

0.98 

1.09 

2.52 

2.30 

2.16 

3C 279 

2.05 

2.06 

2.71 

2.74 

2.73 

PKS 1424+240 

< 0.44 

< 0.50 

< 1.66 

< 1.61 

< 1.56 


TABLE IV: Blazars considered in the Table III For each of them (first column), the deabsorbed value T em 
is reported for different deabsorbing situations. The second column concerns deabsorption according to 
conventional physics, using the EBL model of FRV. The subsequent columns pertain to the photon-ALP 
oscillation scenario with the EBL still described by the model of FRV, and report the values of r em for 
different choices of our benchmark values of the model parameter £ and for Ldom = 4Mpc. The total error 
is the same as for r o b s and reported in the Table III (for more details see text). Sources with question marks 
lack information to perform our analysis and are discarded. For PKS 1424+240 only a lower limit for its z 
exists in the literature: thus, it is neglected in our analysis. 
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Source 

fCF 

L em 

pALP 

1 em 



C = 0.1 

£ = 0.5 

£ = 1 

£ = 5 

3C 66B 

3.00 

3.00 

3.01 

3.01 

3.04 

Mrk 421 

2.05 

2.05 

2.05 

2.07 

2.10 

Mrk 501 

2.54 

2.54 

2.54 

2.57 

2.60 

1ES 2344+514 

2.71 

2.71 

2.72 

2.76 

2.79 

Mrk 180 

3.07 

3.07 

3.08 

3.12 

3.14 

1ES 1959+650 

2.46 

2.46 

2.48 

2.52 

2.55 

1ES 1959+650 

2.32 

2.32 

2.33 

2.38 

2.41 

AP LIB 

2.21 

2.21 

2.23 

2.28 

2.31 

1ES 1727+502 

2.52 

2.52 

2.53 

2.56 

2.58 

PKS 0548-322 

2.44 

2.44 

2.48 

2.56 

2.58 

BL Lacertae 

3.29 

3.29 

3.31 

3.37 

3.39 

PKS 2005-489 

2.77 

2.77 

2.82 

2.89 

2.92 

RGB J0152+017 

2.47 

2.47 

2.53 

2.60 

2.63 

1ES 1741+196 

? 

? 

? 

? 

7 

SHBL J001355.9-185406 

2.81 

2.81 

2.90 

2.99 

3.01 

W Comae 

3.17 

3.17 

3.27 

3.36 

3.38 

1ES 1312-423 

2.17 

2.17 

2.30 

2.38 

2.40 

VER J0521+211 

2.79 

2.79 

2.89 

2.98 

3.00 

PKS 2155-304 

2.81 

2.81 

2.95 

3.03 

3.05 

B3 2247+381 

2.60 

2.60 

2.71 

2.78 

2.80 

RGB J0710+591 

1.89 

1.89 

2.07 

2.15 

2.16 

H 1426+428 

2.85 

2.85 

2.97 

3.06 

3.08 

1ES 1215+303 

2.75 

2.75 

2.92 

3.01 

3.03 

1ES 1215+303 

2.35 

2.35 

2.47 

2.54 

2.56 

RX J1136.5+6737 

? 

? 

? 

? 

7 

1ES 0806+524 

2.70 

2.71 

2.88 

2.98 

3.00 

1ES 0229+200 

0.61 

0.65 

1.41 

1.36 

1.25 

1RXS J101015.9-311909 

2.20 

2.20 

2.41 

2.48 

2.49 

H 2356-309 

2.05 

2.06 

2.30 

2.38 

2.40 

RX J0648.7+1516 

3.45 

3.45 

3.67 

3.75 

3.77 

1ES 1218+304 

1.97 

1.97 

2.27 

2.33 

2.34 

1ES 1101-232 

1.72 

1.73 

2.10 

2.14 

2.13 

1ES 0347-121 

1.87 

1.88 

2.26 

2.29 

2.28 

RBS 0413 

1.88 

1.89 

2.24 

2.31 

2.32 

RBS 0723 

? 

? 

? 

? 

7 

1ES 1011+496 

2.90 

2.90 

3.19 

3.25 

3.26 

MS 1221.8+2452 

? 

? 

? 

7 

7 

PKS 0301-243 

2.93 

2.94 

3.43 

3.48 

3.49 

1ES 0414+009 

1.65 

1.67 

2.25 

2.26 

2.25 

S5 0716+714 

1.60 

1.62 

2.20 

2.22 

2.22 

1ES 0502+675 

? 

? 

? 

7 

7 

PKS 1510-089 

4.02 

4.03 

4.43 

4.47 

4.48 

3C 66A 

1.53 

1.58 

2.40 

2.39 

2.39 

PKS 1222+216 

2.46 

2.48 

2.86 

2.88 

2.89 

1ES 0647+250 

? 

? 

? 

7 

7 

PG 1553+113 

0.98 

1.42 

2.38 

2.22 

2.16 

3C 279 

2.05 

2.13 

2.74 

2.74 

2.73 

PKS 1424+240 

< 0.44 

< 0.70 

< 1.63 

< 1.58 

< 1.56 


TABLE V: Same as Table [TV] but with Ldom = lOMpc. 
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Source 

[cm 

2 s i TeV i ] 

K 

em*' [cm“ 

' 2 s -i TeV i ] 




-^dom 

= 4Mpc 

-^dom - 

= lOMpc 

3C 66B 

1.99 

• lO" 11 

2.13 

•lO" 11 

2.32 

• lO" 11 

Mrk 421 

7.51 

• lO" 10 

8.27 

• lO" 10 

9.18 

• lO" 10 

Mrk 501 

1.80 

• lO" 10 

2.00 

• lO" 10 

2.23 

• lO" 10 

1ES 2344+514 

6.83 

• lO" 11 

7.77 

• lO" 11 

8.80 

• lO" 11 

Mrk 180 

5.82 

•lO” 11 

6.64 

•lO" 11 

7.52 

• 10- 11 

1ES 1959+650 

1.14 

• lO" 10 

1.31 

• lO" 10 

1.49 

1— 1 ■ 

o 

1 

o 

1ES 1959+650 

7.62 

• lO" 11 

8.75 

• lO" 11 

9.94 

• 10- 11 

AP LIB 


? 


? 


7 

1ES 1727+502 

1.29 

• 10- 11 

1.51 

• 10- 11 

1.72 

• 10- 11 

PKS 0548-322 

1.50 

• 10- 11 

1.81 

•10- 11 

2.07 

• 10- 11 

BL Lacertae 

8.83 

• 10" 10 

1.06 

• 10~ 9 

1.21 

• 10” 9 

PKS 2005-489 

4.84 

• 10- 11 

5.86 

•10- 11 

6.70 

• 10- 11 

RGB J0152+017 

2.87 

• 10- 11 

3.53 

• 10- 11 

4.02 

•10- 11 

1ES 1741+196 


? 


7 


7 

SHBL J001355.9-185406 

1.12 

•10- 11 

1.42 

•10- 11 

1.60 

•10- 11 

W Comae 

1.03 

• 10- 10 

1.32 

• 10- 10 

1.47 

I— 1 

o 

1 

o 

1ES 1312-423 

9.03 

• 10- 12 

1.17 

•10- 11 

1.32 

•10- 11 

VER .10521+211 

9.39 

• 10- 11 

1.21 

• 10- 10 

1.33 

• 10- 10 

PKS 2155-304 

2.32 

• 10- 10 

3.04 

• 10- 10 

3.35 

• 10- 10 

B3 2247+381 

2.77 

• 10- 11 

3.58 

• 10- 11 

3.86 

• 10- 11 

RGB J0710+591 

2.49 

• 10- 11 

3.32 

• 10- 11 

3.66 

• 10- 11 

H 1426+428 

2.86 

• 10- 10 

3.76 

• 10- 10 

4.06 

• lO- 10 

1ES 1215+303 

4.38 

• 10- 11 

5.81 

• 10- 11 

6.30 

• 10- 11 

1ES 1215+303 

5.12 

• 10- 11 

6.64 

•10- 11 

6.96 

• 10- 11 

RX J1136.5+6737 


? 


? 


7 

1ES 0806+524 

3.76 

• 10- 11 

5.03 

•10- 11 

5.41 

•10- 11 

1ES 0229+200 

4.04 

• 10- 12 

1.21 

•10- 11 

1.50 

•10- 11 

1RXS J101015.9-311909 

1.51 

•10- 11 

2.03 

•10- 11 

2.17 

•10- 11 

H 2356-309 

2.70 

•10- 11 

3.66 

•10- 11 

3.80 

• 10- 11 

RX J0648.7+1516 

6.83 

• 10- 11 

9.27 

• 10- 11 

9.45 

• 10- 11 

1ES 1218+304 

9.04 

• 10- 11 

1.23 

• 10- 10 

1.24 

• lO- 10 

1ES 1101-232 

4.35 

• 10" 11 

6.26 

• 10- 11 

6.38 

• 10- 11 

1ES 0347-121 

4.37 

• 10- 11 

6.35 

• 10- 11 

6.45 

• 10- 11 

RBS 0413 

3.16 

• 10- 11 

4.41 

• 10- 11 

4.52 

• 10- 11 

RBS 0723 


? 


7 


7 

1ES 1011+496 

1.43 

• 10- 10 

1.92 

• lO- 10 

1.89 

• lO- 10 

MS 1221.8+2452 


? 


7 


7 

PKS 0301-243 

3.95 

•10- 11 

5.47 

•10- 11 

5.32 

•10- 11 

1ES 0414+009 

2.62 

•10- 11 

3.51 

•10- 11 

3.36 

•10- 11 

S5 0716+714 

6.96 

• 10- 10 

9.06 

• lO- 10 

8.63 

• lO- 10 

1ES 0502+675 


? 


7 


7 

PKS 1510-089 

7.75 

•10- 11 

9.87 

•10- 11 

9.37 

•10- 11 

3C 66A 

2.98 

• 10- 10 

4.01 

• lO- 10 

3.81 

• lO- 10 

PKS 1222+216 

1.80 

• 10~ 9 

2.09 

• 10“ 9 

1.98 

• 10” 9 

1ES 0647+250 


? 


7 


7 

PG 1553+113 

4.73 

• 10- 10 

7.86 

• lO- 10 

6.83 

• lO- 10 

3C 279 

2.30 

• 10” 9 

2.28 

• 10~ 9 

2.19 

• 10” 9 

PKS 1424+240 

> 1.92 • 10" 10 

> 2.15 • 10" 10 

> 2.06- lO" 10 


TABLE VI: Blazars considered in the Table III For each of them (first column), the deabsorbed value AT em 
is reported for different deabsorbing situations. The second column concerns deabsorption according to 
conventional physics, using the EBL model of FRV. The subsequent columns pertain to the photon-ALP 
oscillation scenario with the EBL still described by the model of FRV, and report the values of 7v em for 
different choices of our benchmark values of the model parameter Ldom and for £ = 0.5 - which are the best 
values of the model (for more details see text). Sources with question marks lack information to perform 
our analysis and are discarded. For PKS 1424+240 only a lower limit for its z exists in the literature: thus, 
it is neglected in our analysis. 
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Source 

/7>CJP 

r em,AE 

cm 2 s 


77>ALP 

r em,AE 

cm 2 s 

1 




-^dom 

= 4 Mpc 

^dom = 

= 10 Mpc 

3C 66B 

1.75 

■ 10- 11 

1.87 

• io- 11 

2.03 

• io- 11 

Mrk 421 

4.65 

• 10- 10 

5.12 

t— 1 

o 

1 

o 

5.69 

O 

1 

O 

i—i 

Mrk 501 

5.46 

• 10- 11 

6.06 

• 10- 11 

6.76 

• 10- 11 

1ES 2344+514 

2.80 

• 10- 11 

3.19 

• 10- 11 

3.61 

• 10- 11 

Mrk 180 

2.09 

• 10- 11 

2.38 

• 10- 11 

2.70 

• 10- 11 

1ES 1959+650 

3.88 

•10- 11 

4.45 

• 10- 11 

5.05 

• 10- 11 

1ES 1959+650 

2.74 

•10- 11 

3.14 

■ 10- 11 

3.56 

• 10- 11 

AP LIB 


7 


7 


7 

1ES 1727+502 

1.11 

•10- 11 

1.29 

■ 10- 11 

1.48 

• IO" 11 

PKS 0548-322 

2.33 

• 10- 12 

2.79 

■ 10- 12 

3.12 

• IO" 12 

BL Lacertae 

4.41 

• 10- 10 

5.31 

. 10-1O 

6.06 

o 

1 

o 

t-H 

PKS 2005-489 

5.94 

• 10- 12 

7.14 

• 10- 12 

7.99 

• 10- 12 

RGB J0152+017 

4.24 

• 10- 12 

5.17 

• 10- 12 

5.70 

• 10- 12 

1ES 1741+196 


7 


7 


7 

SHBL J001355.9-185406 

7.34 

• 10- 13 

9.09 

• 10- 13 

9.70 

• 10- 13 

W Comae 

1.26 

• 10- 11 

1.59 

• 10- 11 

1.73 

• 10- 11 

1ES 1312-423 

1.42 

• 10- 12 

1.76 

■ 10- 12 

1.82 

• 10” 12 

VER J0521+211 

1.95 

•10- 11 

2.48 

■ 10- 11 

2.69 

• IO" 11 

PKS 2155-304 

5.38 

•10- 11 

6.93 

■ 10- 11 

7.44 

• IO" 11 

B3 2247+381 

1.10 

•10- 11 

1.43 

■ 10- 11 

1.55 

• 10- 11 

RGB J0710+591 

4.79 

• 10- 12 

5.92 

■ 10- 12 

5.81 

• IO" 12 

H 1426+428 

2.04 

•10" 11 

2.66 

■ 10- 11 

2.82 

• 10- 11 

1ES 1215+303 

4.49 

• 10- 12 

5.78 

• 10- 12 

5.96 

• 10- 12 

1ES 1215+303 

3.92 

• 10- 11 

5.17 

• 10- 11 

5.58 

• 10- 11 

RX J1136.5+6737 


7 


7 


7 

1ES 0806+524 

2.86 

• 10- 12 

3.71 

• 10- 12 

3.78 

• 10- 12 

1ES 0229+200 

7.93 

• 10- 12 

6.06 

• 10- 12 

4.84 

• 10” 12 

1RXS J101015.9-311909 

3.10 

• 10- 12 

3.94 

■ 10- 12 

3.89 

• IO" 12 

H 2356-309 

6.02 

• 10- 12 

7.83 

■ 10- 12 

7.75 

• 10- 12 

RX J0648.7+1516 

9.78 

• 10- 12 

1.32 

■ 10- 11 

1.33 

• IO" 11 

1ES 1218+304 

2.85 

•10- 11 

3.67 

■ 10- 11 

3.56 

• IO" 11 

1ES 1101-232 

1.17 

•10- 11 

1.36 

■ 10- 11 

1.22 

• KT 11 

1ES 0347-121 

9.46 

• 10- 12 

1.10 

• 10- 11 

9.86 

• 10- 12 

RBS 0413 

4.64 

• 10- 12 

5.79 

• 10- 12 

5.38 

• 10- 12 

RBS 0723 


7 


7 


7 

1ES 1011+496 

3.93 

• 10- 11 

5.30 

• 10- 11 

5.24 

• 10- 11 

MS 1221.8+2452 


7 


7 


7 

PKS 0301-243 

3.31 

• 10- 12 

4.10 

• 10- 12 

3.79 

• 10- 12 

1ES 0414+009 

7.70 

• 10- 12 

8.49 

■ 10- 12 

7.76 

• IO" 12 

S5 0716+714 

1.69 

• 10- 10 

1.90 

. 10-10 

1.75 

• IO" 10 

1ES 0502+675 


7 


7 


7 

PKS 1510-089 

2.04 

•10- 11 

2.77 

■ 10- 11 

2.69 

• IO" 11 

3C 66A 

3.62 

• 10- 11 

3.58 

■ 10- 11 

3.29 

• 10- 11 

PKS 1222+216 

9.36 

• 10- 10 

1.28 

• 10” 9 

1.26 

• IO" 9 

1ES 0647+250 


7 


7 


7 

PG 1553+113 

1.50 

• 10- 10 

7.59 

• 10- 11 

7.23 

• 10- 11 

3C 279 

9.19 

• 10- 10 

1.14 

• 10“ 9 

1.11 

• IO" 9 

PKS 1424+240 

> 5.65- 10" 11 

> 4.20- 10" 11 

> 4.05 • IO" 11 


TABLE VII: Blazars considered in the Table III For each of them (first column), the deabsorbed value 
Fem,AE is reported for different deabsorbing situations. The second column concerns deabsorption according 
to conventional physics, using the EBL model of FRV. The subsequent columns pertain to the photon-ALP 
oscillation scenario with the EBL still described by the model of FRV, and report the values of F’em.A.E for 
different choices of our benchmark values of the model parameter L & om and for £ = 0.5 - which are the best 
values of the model (for more details see text). Sources with question marks lack information to perform 
our analysis and are discarded. For PKS 1424+240 only a lower limit for its z exists in the literature: thus, 
it is neglected in our analysis. 












